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We derive a simple theory for high-order harmonic generation due to pumping a noninteracting 
metal with a large amplitude oscillating electric field. The model assumes that the radiated light 
field arises from the acceleration of electrons due to the time-varying current generated by the 
pump, and also assumes that the system has a constant density of photoexcited carriers, hence it 
ignores the dipole excitation between bands (which would create carriers in semiconductors). We 
examine the circumstances under which odd harmonic frequencies would be expected to dominate 
the spectrum of radiated light, and we also apply the model to real materials like ZnO, for which 
high-order harmonic generation has already been demonstrated in experiments. 
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I. INTRODUCTION 

The theory of Bloch oscillations has a long history, dat- 
ing back to the original work by Bloch |l| and Zener [2| . 
It remains one of the most difficult phenomena to observe 
in real materials due to the fast timescale for the oscilla- 
tions, and it has been seen predominantly in semiconduc- 
tor superlattices [3[ and in ultracold atomic systems when 
placed on optical lattices [|| . We argue here that similar 
phenomena can be observed by examining the radiated 
light from solids that are pumped by high intensity laser 
pulses with femtosecond pulse widths [a, @ . 

We describe such laser pulses as spatially-uniform 
time- varying electric fields, using a semiclassical ap- 
proach. In this sense, we ignore all magnetic field effects, 
and essentially violate Maxwell's equations for the elec- 
tromagnetic fields. This approach will be accurate if the 
magnetic field effects on the solid are small, which should 
be true for rapidly oscillating fields, because the electron 
spin will not be able to respond to the magnetic field on 
such short time scales. The large electric field will cause 
the electrons in the material to accelerate, giving rise to 
a time-dependent current. Assuming the current on the 
surface is proportional to the bulk current, we will see 
electromagnetic radiation coming off the material with a 
signal proportional to the time derivative of the current. 
We will use this to directly determine the spectra of the 
radiated light. 

Experimental investigations of high harmonic genera- 
tion in bulk crystals are being pursued, including recent 
work on irradiating ZnO crystals, where many tens of 
harmonics of the mid-infrared driving pulse were mea- 
sured [7|. For most crystallographic orientations, only 
the odd harmonics of the fundamental frequency ap- 
peared. We will see below that the occurence of only 
odd harmonics is to be expected for systems that satisfy 
time-reversal symmetry. 



In the next section, we discuss the general formal- 
ism, and apply it to simple model systems, which can 
be thought of as tight-binding models for the bandstruc- 
ture. In Section III, we use numerics to examine the high 
harmonic generation in simple models, and in Section IV, 
we use numerics to examine the high harmonic generation 
in ZnO and compare the numerical results with those of 
experiment. We conclude in the final section. 



II. FORMALISM 

We begin with a system described by a bandstructure, 
which we denote by e(k). The bandstructure is a periodic 
function of the Brillouin zone. If we use a tight-binding 
description for the bands, then the bandstructure can 
generically be written as 



e(k) 



ik<5 



(1) 



for s-orbitals, where d is a translation vector to a neigh- 
boring site (the sum need not stop at just nearest neigh- 
bors), and —tg is the corresponding hopping integral. 
The band index, if there are multiple bands, has been 
suppressed. For models with higher angular momentum 
orbitals, the tight-binding form for the bandstructure is 
more complicated and requires diagonalization of a ma- 
trix whose size is equal to the number of different types 
of orbitals per unit cell (times two in the presence of 
magnetic or spin-orbit interactions). 

We work in the Hamiltonian gauge, where the elec- 
tric field E(f) is described by a spatially uniform vector 
potential A(i), via 



l<9A(t) 
c at 



(2) 



with c the speed of light. The electric field is then intro- 
duced into the bandstructure with the so-called Peierls 
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substitution which takes e(k) — > e(k — eA(t)/hc), with 
e < the electron charge. This substitution properly de- 
scribes the way in which electrons are accelerated in each 
band, but does not take into account the dipole coupling 
between different types of bands (like between s and p 
bands). Such a direct Zener tunneling term can be in- 
cluded, but will complicate the analysis, and we ignore it 
in the results we describe here, for simplicity. We expect 
that the results one finds if we include such terms will 
be similar, because the main effect of the laser pulse is 
to accelerate the electrons within a given band in metal- 
lic systems, and excitations between bands become more 
important for semiconductors or insulators. But even in 
that case, if we assume most of the photoexcitation oc- 
curs when the field amplitude is the highest, then the 
results of this work will be identical except for a factor 
of 2, which does not affect any of the results we present 
below. 

The Hamiltonian of the system then becomes 



If we have a pulse with a shape given by 



kcr 



where (c ko .) is the electron creation (annihilation) op- 
erator for an electron with momentum k and spin a (here 
we are assuming, for simplicity, a single band model and 
no spin-orbit coupling, so a is a good quantum number) . 
This is a time-dependent Hamiltonian, but it commutes 
with itself at different times, which makes its time evo- 
lution particularly easy to solve. 

The system starts off in equilibrium at early times t — > 
— oo and the electric field pulse is centered around t = 0. 
Hence, the initial distribution of the electrons is given by 
the Fermi-Dirac distribution 



/(e-M) 



1 + exp[/?(e - n)] 



(4) 



in the distant past. Here, /i is the chemical potential, and 
(3 is the inverse temperature. We work with an isolated 
system, so the time dynamics can be solved directly in 
the Heisenberg picture. The current satisfies 



krr 



eA(t) 

he 



(5) 



which simply sums the field-dependent velocity (Vke/ft), 
weighted by the initial distribution of the electrons and 
multiplied by the electric charge. This formula also fol- 
lows from a Kadanoff-Baym-Keldysh formalism Q where 
one can exactly solve for the Green's functions, evaluate 
them at equal time and sum the electron occupancy at k 
weighted by the field-dependent velocity. 

If we assume that the radiated light arises from the 
acceleration of the electrons, then the light spectrum is 
proportional to 



dte lult —3{t) 



dte iwt 3(t) 



(6) 



A(t) = — Aq cos(cjq£ 



exp^ 



(7) 



then we are guaranteed that the field has no dc compo- 
nent (required of optical pulses) and if to is somewhat 
larger than the period of the oscillations 2tt/uiq, then the 
pumped electric field has a reasonably well defined fun- 
damental frequency (4> is an adjustable phase shift of the 
oscillations relative to the gaussian envelope) . While we 
choose a gaussian form for the envelope function, any nor- 
malized slowly varying function which vanishes rapidly 
enough as \t\ — > oo can be used for the envelope, which 
we denote below by g(t). By calculating the spectrum 
of the radiated light from the Fourier transform of the 
current as a function of time, we can examine the high 
harmonic generation in the system. 

Before jumping into numerics, we first discuss some 
general principles and examine symmetry considerations. 
If the lattice is inversion symmetric, then the parity op- 
eration immediately tells us that the energy is an even 
function of the wavevector 



e(k) = e(-k). 



(8) 



If the lattice lacks inversion symmetry, this result still 
holds if there is no spin-orbit coupling and no magnetic 
field due to magnetic order. This is because the Hamil- 
tonian is time-reversal symmetric, so there is a Kramers 
doublet, 



e t (k) = e 4 (-k). 



(9) 



Because the system is also degenerate with respect to 
spin [since the operator Si commutes with H(t)], we 
find e-|-(k) = e^(k), which then implies Eq. ((5J). Finally, 
we also examine what happens in the presence of spin- 
orbit coupling of atomic origin. In this case, since the 
spin-orbit potential has the periodicity of the lattice, the 
eigenstates remain in the Bloch form, but since each com- 
ponent of spin no longer commutes with the Hamiltonian, 
the spin becomes entangled with the wavevector. We con- 
tinue to have a Kramers doublet, where the spin label in 
Eq. © is replaced by a pseudospin. Because k — > k 
and f - ^4- under a time-reversal operation, the two de- 
generate states will have opposite velocities. Hence, in 
equilibrium (A = 0), the total current vanishes because 
the contributions from the states with opposite momen- 
tum and opposite pseudospin cancel in the summation in 
Eq. ©. 

Now consider the case where Eq. © holds (inversion 
symmetry or time-reversal invariance and no spin terms 
in the Hamiltonian). Then, we can replace the exponen- 
tial in Eq. ([1} by a cosine, so we have 



e(k) =-^> cos(k-<5). 



(10) 
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The field-dependent velocity then becomes 



v k (t) 



-y 



St. 



s sva 



eA(t) 



he 



(11) 



Using the trigonometric summation formula, we can ex- 
pand the sine function into the difference of products of 
sines and cosines with arguments of k • 8 and A • 8. Sub- 
stitute this into the equation for the current [Eq. ([5])] 
and recall that the bandstructure is an even function of 
k. This implies that the only term that survives in the 
summation is the term proportional to cosk- 8, since the 
sine term is an odd function of the wavevector. This then 
yields 

J ( < ) = \ E E cos ( k • 5 ) sin f^T^) A*)- 

k s \ c J 

(12) 

After performing the summation over momentum, the 
current will be expressed as a constant amplitude for each 
8, 3(8), multiplied by the sine function which contains 
all of the time dependence, and then summed over the 
neighbors. We manipulate the sine function as follows: 

'eAftl • 8 

" 1 = -sin[E • 8cos(uj t + (j))g(t)} (13) 



sm 



he 



= -ImJ J [B -8g(t)} 



(14) 
(15) 



2 J «I E o ' Sg(t)] cos[n(uj t + (/>)] 



71=1 



2^(-l)»J 2n+1 [Eo-Js(t)] 

n=l 

x cos[(2n + l)(u t + (/>)]. 



(16) 



Here, we use the general formula g(t) for the slowly vary- 
ing envelope function, and the symbol J n denotes the 
nth Bessel function. The Fourier transform of this will 
correspond to sharp peaks around the odd harmonics 
(2n + l)wo, with the width of the peak determined by 
how rapidly the envelope function g(t) decays in time. 
The amplitude of the peak can decrease rather rapidly, 
as the Bessel functions decay for large index, with a fixed 
argument. Hence, one expects to generically see the odd 
harmonics in the high harmonic generation if the system 
obeys the relationship in Eq. ((5J). It is straightforward 
to generalize this argument to non-inversion-symmetric 
systems with spin-orbit coupling if the total current is 
simply the sum over contributions from the two bands 
with different pscudospin. 



III. NUMERICAL RESULTS ON MODEL 
SYSTEMS 

We will examine the generic current response for hy- 
percubic lattices in d-dimensions, where the field Eq has 
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FIG. 1: Top panel: normalized time-dependent electric field 
with ojo = 2, cj> — it/2, and to = 10 in units of hwo/ea. Bottom 
panel: current when Eq = 10 in units of B. 



equal components along each coordinate direction, or has 
no component along a coordinate direction [hence we ex- 
amine (1,0,0), (1,1,0), or (1,1,1) directions for a simple 
cubic lattice in d = 3] . Using the fact that all coordinate 
directions are equivalent, so ^ k cos kiaf(e\ i — //) is inde- 
pendent of the coordinate direction i (with a the lattice 
constant), then gives that the components of the current 
satisfy 



Ji(t) 



^-^sinA 4 (t) / deep(e)f(e 



M), (17) 



where p(e) = ^ k <5(ek — e) is the noninteracting density of 
states of the d-dimensional hypercubic lattice. This form 
of the current involves a temperature-dependent ampli- 
tude term B — ea J deep(e)f(e — p)/2dh and the time- 
dependent term 



Jj(i) / he 

—— = sm E 0i — cos(w i 
B V ea 



1 



-(t/t ) 



(18) 



Here we define Eq = A^hct^/eci which is the dimension- 
less amplitude of the electric field (the units of the field 
are fuvo/ea). The radiated light is then proportional to 
\u}Ji(ui)\ 2 , which we find by taking the Fourier transform 
of Eq. (|18[) . Note that this specific form of the radiated 
harmonics is independent of both the direction of the field 
and the dimensionality of the lattice for noninteracting 
electrons on hypercubic lattices. 

We begin by showing the normalized electric field in 
the top panel of Fig. Q] We choose ujq = 2, </> = tt/2, and 
to = 10. This produces an even (in time) field profile 
that has on the order of 10-12 oscillations clearly visible, 
allowing for the frequency of the pump pulse to be rea- 
sonably well defined. The current as a function of time 
(without the temperature and dimension-dependent am- 
plitude factor B) is shown in the bottom panel. One can 
clearly see that as the field amplitude rises, the periods of 
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FIG. 2: Spectrum of radiated light from an electric field pulse with ujo = 2, cj> — n/2, to = 10, and Eo taking on a range of 
different values: (a) E = 2, (b) E = 5, (c) E = 10, (d) E = 20, (e) E = 40, and (f) E = 80. Note how the number 
of harmonics increases as the field amplitude increases, and how the signal appears to become noisier at small frequencies for 
higher field amplitudes. 



the oscillations in the current get shorter, and that there 
is a "clipping" like effect in the oscillations when the 
field amplitude is large enough because the sine function 
is bounded between plus and minus one. The current has 
significant structure in it as a function of time and hence 
will give rise to a complicated Fourier transform. But, by 
just looking at the results in Fig. [TJ we can reason that 
we expect to see higher frequencies in the Fourier trans- 
form when the field amplitude is higher. Indeed, this is 
what we next show. 

In Fig. [21 we plot the logarithm of the radiated light 
spectrum in arbitrary units, as a function of the fre- 
quency of the light. The field profile is identical to that 
found in Fig. [TJ and the different panels show the ef- 
fect of increasing the amplitude Eq of the pump pulse. 
The plot is a semilogarithmic plot, because the higher 
harmonics are often suppressed by orders of magnitude 
from the fundamental before we hit a "noise floor" in 
the data. Even though we have a completely well de- 
fined integrand to integrate over, and we use an adaptive 
integration routine available from quadpack for handling 
such Fourier integrals, it still is challenging to obtain high 
accuracy for these integrations. The main source of er- 
ror comes from balancing the integration range versus 
the order of the polynomials used in the subdivided in- 
tegration ranges and the absolute accuracy requested in 
the subroutine. We believe that most of our results with 
an amplitude larger than 10 -10 are completely accurate, 
but we do not know whether there might be additional 
structure in the noise region. Note that this is purely 
a numerical issue, and will have little bearing on experi- 
ment, because once the higher harmonics are too small in 
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FIG. 3: Fine structure for low frequencies when wo = 2, <f> = 
7r/2, to — 10, and Eo = 20 (magnification of panel [2] (d) 
above). Note how the extra structure becomes less frequent 
for higher frequencies. 



amplitude, they cannot be measured. In general, we see 
on the order of a few up to maybe 20 higher harmonics 
appearing clearly in the figures. 

One can see two main trends occurring as the ampli- 
tude is increased: (i) first, the start of the noise floor 
moves to higher frequencies as the field amplitude is in- 
creased, as we expected, and hence we can see more 
higher harmonics as the system is driven harder, and 
(ii) we see first the development of a flat plateau of high 
amplitude harmonics occur at lower frequencies, which 
increases in size as the field amplitude increases. But, 
we also see the development of additional structure that 



looks like noise for these high field-amplitude curves. 
To investigate this further, we have enlarged the low- 
frequency region of panel [2] (d) in Fig. [3] for the case 
with Eq = 20. One can see that this fine structure is 
not noise, as it appears in the figures with the wider fre- 
quency range, but is instead a fine structure of additional 
oscillations in the harmonic peaks that show more struc- 
ture the lower the harmonic frequency is. The structure 
is reminiscent of a beating effect, but we have not been 
able to completely track down its origin. 

In Fig. 21 we investigate the effect of the phase of 
the field on the high harmonic generation. The effect 
is rather mild, as can be seen in the figure, where the 
different colors correspond to different phases, but the 
evolution is definitely not monotonic in the phase. Since 
detailed field profiles and phases of the electric field are 
difficult to determine for a given pump pulse, this weak 
phase dependence makes the experimental results much 
more robust, since the most important aspects of the field 
are the pulse width, oscillation frequency, and amplitude. 
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FIG. 4: (Color online) Effect of the phase (j> on the fine struc- 
ture. The different colors correspond to different phases: 
black 0/2); red (3tt/8); green 0/4); blue 0/8); and n S ht 
blue (0). Note how the dependence is nonmonotonic in the 
phase, but overall is quite mild. 



IV. NUMERICAL RESULTS FOR ZINC OXIDE 

We now apply this model to ZnO using the band 
structure calculated within density functional theory 
(DFT). The wurtzite crystal structure of ZnO is uniaxial 
and non-centrosymmetric, and spin-orbit coupling only 
weakly affects the band structure (the splitting between 
bands at the top of the valence band is only a few meV) . 
Thus based on the symmetry arguments presented in Sec. 
II, we expect to see only odd harmonics in the current 
response. 

The DFT calculations were performed usingthe Vi- 




enna Ab Initio Simulation Package (VASP) with 
the electron-ion interaction described using projector 
augmented wave potentials [l2j], and the exchange- 



FIG. 5: Band structure of wurzite ZnO calculated within 
GGA-DFT. The zero of energy is set to the top of the valence 
band. The arrows indicate the chemical potentials needed for 
a 10% electron filling of the conduction band and a 10% hole 
filling of the valence bands. 



correlation interaction handled with the generalized gra- 
dient approximation (GGA) [l3| . Spin-orbit effects were 
not included. The structural parameters were optimized, 
and the Kohn-Sham eigenvalues were tabulated on uni- 
form grids in the Brillouin zone. Band velocities were 
calculated using centered differences. We have found 
that the noise floor in the calculated spectrum converges 
slowly with respect to the Brillouin-zone sampling. With 
a grid of 64 x 64 x 32 k-points, the noise floor in \u>J(oj)\ 2 
lies about 7 orders of magnitude below the maximum, 
which is sufficient to distinguish on the order of ten or 
more harmonics, depending on the magnitude of the field. 

The calculated band structure is shown in Fig. [5] The 
highest lying valence bands have mostly O 2p character, 
while the conduction band is primarily of Zn 4s charac- 
ter. The GGA severely underestimates the ZnO band 
gap, yielding a value of about 0.8 eV as compared to the 
measured gap of 3.4 eV. Quasiparticle calculations based 
on the GW approximation have found that the main cor- 
rection to the DFT bands for ZnO is an upward shift by 
about 2 eV of the two lowest conduction bands [14| . The 
overall shape and width of the O 2p valence bands and of 
the Zn 4s conduction bands do not change significantly. 
Since we address here only the contribution to high-order 
harmonic generation from field-induced acceleration of 
electrons within bands, and ignore contributions from 
field-induced dipole transitions between bands, the GGA 
bands provide a reasonable starting point for a band-by- 
band analysis. 

Because we are not modeling the photoexcitation pro- 
cess, we need to artificially shift the chemical potential 
to create carriers in this semiconducting system. In the 
results presented below, the response of electrons in the 
conduction band is calculated by setting the chemical po- 
tential to the equilibrium value corresponding to a 10% 
occupation of the lowest conduction band. Given the 
large width of the Zn 4s band, this corresponds to a 
chemical potential about 2 eV above the conduction band 
minimum, as marked in Fig. [5] The O 2p bands have 
larger multiplicity and are much narrower than the con- 
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FIG. 6: Spectrum radiated by carriers in the ZnO conduction 
band when driven by a pulse with the electric field polarized 
along [0001]. The maximum electric field strength is (a) 0.145 
V/A, (b) 0.29 V/A, and (c) 0.58 V/A. 
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FIG. 7: Effects of field direction and band dispersion on the 
emitted spectrum. Panels (a) and (b) show emission by elec- 
trons in the ZnO conduction band when driven by out-of- 
plane and in-plane electric fields, respectively. Similarly, pan- 
els (c) and (d) show emission by holes in the ZnO valence band 
when driven by out-of-plane and in-plane fields, respectively. 



duction band, so a 10% hole occupancy of valence bands 
requires a much smaller shift in the chemical potential 
with respect to the band edge. Without a mechanism 
for interband transitions, the response from holes is es- 
timated by simply summing the current from all the va- 
lence bands with partial occupation. 

A pulse of the form given by Eq. (8) is applied, with 
phase = 0, width t = 40 fs, frequency uj = 20/i , and 
maximum electric field strengths on the order of 0.1 to 1 
V/A. These parameters yield pump pulses similar to the 
~ 9 cycle, mid-infrared laser pulses used in the recent ex- 
perimental study of high-harmonic generation in ZnO 0] ■ 
Figure [5] shows the spectrum of light emitted by electrons 
in the conduction band when driven by an electric field 
polarized along the [0001] direction. The maximum field 
strength of the pulse doubles between successive panels, 
showing the effect of field amplitude on the spectrum. 
The same general trends that were observed in the cal- 
culation based on tight-binding bands are evident here. 
In particular, when the the field amplitude is increased, 
a plateau develops in the low- frequency regime, and the 
low-order harmonic peaks acquire more and more fine 
structure. In the experimental setup in Ref. 0] the de- 
tection cut-off at low-energies was set approximately to 
the band gap, so the reported spectra start from about 
the 9th order harmonic and range up to about the 25th 
harmonic. Within this range, the calculated decay of the 
intensity with increasing harmonic order agrees roughly 
with the measured spectra. It would be interesting to try 
to extend the measurements to lower-order harmonics to 
see if the fine structure that appears in our calculations 
is observed. 

Given the hexagonal symmetry of the wurtzite crystal 
structure, it is expected that the response of electrons in 
ZnO will depend on whether the electric field is paral- 



lel or perpendicular to the c axis. In addition, since the 
dispersion of the conduction and valence bands are quite 
different, it is interesting to compare the emission spectra 
for electrons and holes. Figure [7] shows the spectra gen- 
erated by electrons in the conduction band and by holes 
in the valence band when driven by fields in the [0001] 
and [llOO] directions. The parameters for the pump pulse 
(with the exception of the field direction) are the same 
as those used for Fig. [B^. When the driving field lies 
in the basal plane (Fig. [TJd and d), the low- frequency 
plateau region is shortened compared to when the field 
is perpendicular to the plane (Fig. [7^, and c). The in- 
plane field produces a stronger signal for the lower har- 
monics, but the magnitude of the harmonic peaks decays 
more quickly with harmonic order. Further, the low- 
order peaks are narrower and exhibit less structure than 
when the field is polarized perpendicular to the plane. 
Despite significant differences in e(k) for the electrons 
and holes, the spectral characteristics of their radiation 
are very similar. There is an overall reduction in signal 
intensity from the valence band as compared to the con- 
duction band, which is consistent with the holes having 
a signficantly larger band mass than the electrons. Inter- 
estingly, the fine structure in the lower harmonics seems 
to depend much more on the direction of the electric field 
than on whether the carriers are in the valence or con- 
duction band. At higher frequencies, the detailed shape 
of the peaks exhibits stronger dependence on the band 
structure. Varying the filling factors for the conduction 
and valence bands changes some of the detailed features 
in the spectra as well as the overall intensity of emitted 
light, but the general trends discussed here seem to be 
generic and robust. 

A remaining puzzle is that while most crystallographic 
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orientations produced only odd-order harmonics in the 
experiments, both even and odd harmonics were observed 
when the polarization of the incident electric field was 
parallel (or close to parallel) to the c axis of the ZnO 
crystal [7j. Our model based on intraband dynamics of 
carriers in the bulk is not able to account for the ap- 
pearance of even harmonics under these circumstances. 
To more fully understand the experimental results for 
ZnO (or other semiconductors), the coupling between 
bands @ and other effects such as the broken symme- 
try at the surface need to be further explored. 



on the Kadanoff-Baym-Keldysh contour [j| [T(|, which 
we have carried out an initial study for elsewhere 
It will be interesting to examine how this behavior sur- 
vives as correlations increase, into materials like Mott 
insulators. In addition, it will be important to under- 
stand the balance between the more universal behavior, 
which is captured by model system calculations, and the 
materials specific features, which require more materials 
specific calculations based on density functional theory 
plus extensions to include correlations in the nonequilib- 
rium domain. The future for these types of studies, both 
theoretical and experimental is clearly a bright one! 



V. CONCLUSIONS 



In this work, we have shown some simple theories for 
high harmonic generation in metals (or photodoped semi- 
conductors) which capture many of the experimental re- 
sults already seen in this process. While it has long been 
thought that Bloch oscillations occur at too fast a time 
scale to be observed experimentally in metals, this par- 
ticular experiment is the exception, which allows one to 
indirectly observe the oscillations from the light gener- 
ated by the accelerating electrons. Of course, in any real 
system, one needs to include interactions to determine 
how feasible the observation of these effects are. The 
most important ones will be interactions with defects 
(disorder), with phonons, and with other electrons as the 
correlations in the material grow. In order to do this, a 
more complete many-body physics formulation is needed 
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